function M = getM(n_seg, n_order, ts)
    M = [];
    n_poly_perseg = n_order + 1;
    d_order = n_poly_perseg / 2;
    for k = 1:n_seg
        %#####################################################
        % STEP 1.1: calculate M_k of the k-th segment 
        %
        %
        %
        %
        M_k = zeros(n_poly_perseg, n_poly_perseg);
        for i = 1 : n_poly_perseg
            if i <= d_order
                t = 0;
                cur_d=i-1;
            else
                t = ts(k);
                cur_d=i-1-d_order;
            end
            for j=cur_d:n_order
               M_k(i,j+1)=factorial(j)/factorial(j-cur_d)*t^(j-cur_d);
            end
        end
        
        M = blkdiag(M, M_k);
    end
end